import sklearn
import numpy as np
from scipy import ndimage

from matplotlib import pyplot as plt
%matplotlib inline
import matplotlib.patches as mpatches

import skimage
from skimage import io
from skimage.filter import threshold_otsu

from skimage.feature import peak_local_max
from skimage.morphology import watershed

img1 = skimage.io.imread("img/cyano1.jpg")

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (9,9) )
sub1 = plt.imshow(img1)

arrow_args = dict(arrowstyle="->")

sub1.axes.annotate('Heterocyst'
                   ,xy= (1600,700)
                   ,xytext = (700,300)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 14
                   )
sub1.axes.annotate('Vegetative'
                   ,xy = (1080,960)
                   ,xytext = (400,600)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 14
                   )

def boundbox1(tryit):
#     np.nonzero(np.apply_over_axes( np.argmax, thelabel,[1]))[0].min()
    minX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].min()
    maxX = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[1]))[0].max()
    minY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].min()
    maxY = np.nonzero(np.apply_over_axes(np.argmax,(tryit),[0]))[1].max()
    return [minX,maxX,minY,maxY]

def buildrect(arr):
    
    minY,maxY,minX,maxX = arr
    
    rect1 = mpatches.Rectangle(xy=(minX,minY)
                           ,width = maxX - minX
                           ,height = maxY - minY
                           ,fill = False
                           ,edgecolor = 'red'
                           ,linewidth = 2)
    return rect1

def img_to_labels1(inp_img):
    
    thresh = threshold_otsu(inp_img)
    otsu0 = inp_img < thresh
    dist0 = ndimage.distance_transform_edt(otsu0)
    j=34
    peaks0 = peak_local_max(dist0 ,indices = False
                            ,footprint=np.ones((j,j))) 
    print peaks0.cumsum()
    marker0 = ndimage.label(peaks0)[0]
    labels0 = watershed(-dist0, marker0, mask=otsu0)
    
#     plt.imshow(labels0)
#     plt.show()
    return labels0

def make_windows(inp_labels, inp_img):
    bb0 = []
    for i in range(1,inp_labels.max() + 1):
        label_i = (inp_labels == i)
        bb0.append(boundbox1(label_i))

    rect_arr = [buildrect(bb0[i]) for i in np.arange(bb0.__len__())]

    fig = None
    s0 = None
    fig, (s0) = plt.subplots(nrows = 1, ncols = 1, figsize = (12,12))
    s0.imshow(inp_img, cmap='gray')
    for rect in rect_arr:
        s0.add_patch(rect)
    fig.show()
    return fig

def algo1(inp_img):
    labels = img_to_labels1(inp_img)
    out_fig = make_windows(labels,inp_img)
    return out_fig

img1_gray = skimage.color.rgb2gray(img1)

plt.imshow(img1_gray, cmap='gray')

cross = img1[1200:2000,0:500]

cross = skimage.color.rgb2gray(cross)

algo1(cross)

inp_img = cross
thresh = threshold_otsu(inp_img)
otsu0 = inp_img < thresh
dist0 = ndimage.distance_transform_edt(otsu0)
j=34
peaks0 = peak_local_max(dist0 ,indices = False
                        ,footprint=np.ones((j,j))) 
print peaks0.cumsum()

img = peaks0
peaks_array =[]

peak_counter = 0

for y in range(0,img.shape[0]):
    for x in range(0,img.shape[1]):
        
        if (img[y,x] == True):

            peak_counter = peak_counter + 1
            peaks_array.append([peak_counter,y,x])
            

img = cross
offset = 40

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')

arrow_args = dict(arrowstyle="->")

for i in peaks_array:
    
    sub1.axes.annotate('Peak '+ str(i[0] - 1)
                   ,xy= (i[2],i[1])
                   ,xytext = (i[2]+offset,i[1]+offset)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 12
                   )

X = np.array(peaks_array)[:,(1,2)]

from sklearn.neighbors import kneighbors_graph

weighted_graph = kneighbors_graph(X, 5, mode='distance')

def strand_v1(dist, pres, past):
    """
    args:
        dist: sklearn distance kneighbors_graph
        pres: present label
        past: past label
    """
    m, n = dist.shape
    n1nbr = None
    for i in range(n):
        if i in past or dist[pres, i]==0:
            continue
        if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
            n1nbr=i
#     print pres, past, n1nbr
    if len(past) == 0 or n1nbr is not None:   
        return strand_v1(dist, n1nbr, past+[pres])
        
    return past + [pres]

def strand_v2(dist, pres, past):
    """
    Uses a mean std deviation metric to ensure neighbors are cells reasonably close
    to one another, based on 
    args:
        dist: sklearn distance kneighbors_graph
        pres: present label
        past: past label
    """
    m, n = dist.shape
    n1nbr = None
    for i in range(n):
        if i in past or dist[pres, i]==0:
            continue
        if n1nbr is None or dist[pres,i]<=dist[pres,n1nbr]:
            n1nbr=i
#     print pres, past, n1nbr
    if len(past) == 0 or n1nbr is not None: 
        rawr = dist[]
        if n1nbr <  and np.std(past+[pres])==0:
            return strand_v2(dist, n1nbr, past+[pres])
        
    return past + [pres]

rawr = np.random.rand(10)*500

print rawr

for i in range(1, len(rawr)):
    print mean_std(rawr[:i])

strand_v2(weighted_graph, 3, [])

strand_v2(weighted_graph, 19, [])

strand_v1(weighted_graph, 15, [])

strand_v1(weighted_graph, 5, [])

hand_labeled = [19, 13, 12, 10, 9, 7, 8]
# labels=tuple(strand_v1(weighted_graph, 19, []))
labels = hand_labeled
print labels

X[labels,:]

coefficients = np.polyfit(X[labels,1], X[labels,0], len(labels)-1)
print coefficients

polynomial = np.poly1d(coefficients)

ys = polynomial(X[labels,1])

fig,ax=plt.subplots(figsize=(15,15))
ax.imshow(img, cmap='gray')
ax.plot(X[labels,1], ys, 'r', lw=3)
ax.set_ybound(img.shape[0],0)
ax.set_xbound(img.shape[1],0)

img = cross
offset = 40

bplt, (sub1) = plt.subplots(nrows=1,ncols=1,figsize = (25,25) )
sub1 = plt.imshow(img, cmap='gray')

arrow_args = dict(arrowstyle="->")

for i in peaks_array:
    
    sub1.axes.annotate('Peak '+ str(i[0] - 1)
                   ,xy= (i[2],i[1])
                   ,xytext = (i[2]+offset,i[1]+offset)
                   ,xycoords = 'data'
                   ,textcoords = 'data' #'offset points'
                   ,arrowprops = arrow_args
                   ,fontsize = 12
                   )

# east_strand = strand_v1(weighted_graph, 19, [])
east_strand = hand_labeled

north_bound = east_strand + [11, 6]
print north_bound

west_bound = east_strand + [11, 14]
print west_bound

south_bound = east_strand + [11, 17]
print south_bound

def show_path_and_der(coords, labels, do_plot=True):
    labels = tuple(labels)
    coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
    polynomial = np.poly1d(coefficients)
    ys = polynomial(coords[labels,1])
    dcoefficients = np.polyder(coefficients)
    print dcoefficients
    dpolynomial = np.poly1d(dcoefficients)
    dys = dpolynomial(coords[labels,1])
#     d2coefficients =np.polyder(coefficients, 2)
#     print d2coefficients
#     d2polynomial = np.poly1d(d2coefficients)
#     d2ys = d2polynomial(coords[labels,1])
    
    if do_plot:
        fig,ax=plt.subplots(figsize=(10,10))
        ax.imshow(img, cmap='gray')
        ax.plot(coords[labels,1], ys, 'r', lw=3)
        ax.plot(coords[labels,1], dys, 'g', lw=3)
#         ax.plot(coords[labels,1], d2ys, 'b', lw=3) 
    #     ax.set_ybound(img.shape[0],0)
    #     ax.set_xbound(img.shape[1],0)

show_path_and_der(X, north_bound[:-1])
show_path_and_der(X, north_bound)

show_path_and_der(X, west_bound[:-1])
show_path_and_der(X, west_bound)

show_path_and_der(X, south_bound[:-1])
show_path_and_der(X, south_bound)

def cost_of_path(coords, labels):
    labels = tuple(labels)
    def get_der_coef(coords, labels):
        coefficients = np.polyfit(coords[labels,1], coords[labels,0], len(labels)-1)
        dcoefficients = np.polyder(coefficients)
        return dcoefficients
    initial = get_der_coef(coords, labels[:-1])
    decision = get_der_coef(coords, labels)[1:]
    import scipy
    return scipy.spatial.distance.euclidean(initial, decision)

cost_of_path(X, north_bound)

cost_of_path(X, west_bound)

cost_of_path(X, south_bound)